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Abstract 


We report Markov chain Monte Carlo fits of the thermophysical model of Wright (20071 to the fluxes of 10 
asteroids which have been observed by both WISE and NEOWISE. This model is especially useful when 
one has observations of an asteroid at multiple epochs, as it takes advantage of the views of different local 
times and latitudes to determine the spin axis and the thermal parameter. Many of the asteroids NEOWISE 
observes will have already been imaged by WISE, so this proof of concept shows there is an opportunity to 
use a rotating cratered thermophysical model to determine surface thermal properties of a large number of 
asteroids. 
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1. Introduction 


Thermophysical asteroid models have been in use for decades, and they have gradually improved in their 
completeness. The effects of thermal inertia in causing a time delay from local noon in the temperature max¬ 
imum have long been taken into account ( {Peterson 1976). The peaking of emission near zero observational 
phase angle (‘beaming’) was accounted for later in the Standard Thermal Model (STM) by calculating the 
emission at zero phase and then applying a linear correction factor for other phases. A beaming correction 
factor was then applied to account for the fact that beaming reduces reradiated energy; the STM is from 


Lebofsky et al. (1986). The beaming correction had also been used in Jones & Morrison (1974), and Morri¬ 
son & Lebofsky ( [1979 ). Harris (1998) improved the effectiveness of the STM with his Near Earth Asteroid 
Thermal Model (NEATM) by letting this correction factor be a free parameter, and adjusting it to match 
the observed color temperature when multiple thermally-dominated wavelengths are available. However, 
both of these are empirical models that account for a variety of phenomena and parameters using only the 
beaming parameter. 


Hansen (1977) did not use a beaming model and instead considered an asteroid as covered in craters, so 


that at non-zero phase angle there is increased shadowing of the visible portion of the asteroid over what 
would be observed from a smooth surface. This dampening of flux at non-zero phase can be interpreted as a 


peak at zero phase. Spencer (1990) adds to this model the consideration of light reflecting off different parts 


of craters, as well as an iterative numerical process to model heat conduction. Lagerros (1996) combines 


the effects of both thermal inertia and cratering in his model, and calculates a correction factor based on 
a comparison between a smooth surface and one with craters, though a more detailed discussion as to the 


effects of different sorts of surface roughness is given in Lagerros (1998). Delbo’ et al. (2007) includes surface 


roughness by considering the mean slope of the surface of the asteroid, rather than assuming any sort of 


crater geometry. Hanus (2015) uses optical photometric data to investigate shape models of the asteroid 


before applying a thermophysical model using infrared data. 
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Wright (2007) takes the surface cratering into account explicitly in his Spherical, Cratered, Rotating, 


Energy-conserving Asteroid Model (SCREAM, name assigned here for ease of reference) by including the 
local effects of this geometry in the power balance calculations of the temperature distribution over the 
surface of the asteroid. As a result energy is entirely conserved, and the model can include the effects of the 
reflection of solar light and the absorption of blackbody radiation caused by the mutual visibility of different 
parts of a crater, in addition to considering vertical heat conduction. Other asteroid thermophysical models 


that take all of these into account include those of Muller (2007), Rozitis & Green (20111, and Leyrat et al. 

The Wide-held Infrared Survey Explorer (WISE) mission has provided a veritable treasure trove of 
information on the infrared sky (Wright et al.[ 2010). This includes asteroids, of which over 160,000 have 


now been observed. The NEOWISE project allowed individual exposures from WISE to be publicly archived 
and searched for moving objects, to enable the discovery of new asteroids and comets (Mainzer et al. 2011a). 

The NEATM makes it possible to quickly perform thermal modeling of asteroids and has already been 
used on WISE data (e.g. Mainzer et al. 2011c). While the SCREAM is much more computationally 


intensive, it has the potential to allow additional parameters beyond diameter, albedo, and beaming to be 
determined, such as spin axis and thermal inertia. Parameters such as thermal inertia and spin axis can be 
more narrowly constrained when observations of an asteroid are available at multiple epochs. In cases when 
multiple viewing geometries are available, the NEATM can converge to different beaming factors at each 
epoch (though this can also result just from asphericity or different viewing geometries). Multiple epochs 
of observation are very advantageous in the SCREAM, as the differing phase angle gives views of different 
local times and/or latitudes of the asteroid, which allows one to characterize the asteroid spin axis in order 
to explain the phase-varying flux. 

With the recent reactivation of the WISE telescope for the restarted NEOWISE mission, many asteroids 
are now being reobserved at different phase angles (Mainzer et al. 2014). This new mission thus gives us an 
opportunity to characterize these asteroids using the SCREAM, with which we can jointly fit all the data 
to explain the phase-varying flux. As a proof of concept, we here report Markov chain Monte Carlo fits of 
the SCREAM to 10 asteroids which have already been reimaged by the NEOWISE mission. 


2. Data 

Candidates for analysis were found by querying the Minor Planet Centeij^ (MPC) for all WISE and 
NEOWISE observations of asteroids, and then searching through the output for asteroids which were seen 
by both. Then the Infrared Science Archive’^ moving object search feature was used to find flux data 
for each of the asteroids. After throwing away temporal outliers (> 1 day from other observations), the 
data were binned into time series from different observational epochs, and then the interquartile mean (or 
mid-average) of each epoch was taken as the new data point. Since our asteroids all have prograde orbits 
with periods > 1 yr, no meaningful intra-bin trends were seen in the observational epochs, which were of 
length < 10 days. A new uncertainty for each data point was calculated as: 
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where there are N sorted observations being mid-averaged with uncertainty aj, and their mid-average flux 
value is fi. The first term is the standard combination of independent uncertainties applied to the second 
and third quartiles of the data, but the correction factor C ~ 0.77 accounts for the extra information from 


^ htt p: / / WWW .minorplanetcenter.net/ 
^https://irsa.ipac.caltech.edu/frontpage/ 
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the data points we discarded and was found via Gaussian error modeling The second term is the equivalent 
of 0.1 magnitudes, and was added to account for the magnitude of the approximations made in our model 
which are detailed in Section especially the discretization of the craters. 

A summary of the data used can be found in Table 

When many high-accuracy observations of an asteroid over its rotational period are available, one may 
use a technique called ‘lightcurve inversion’ to deduce both the shape and rotational characteristics of the 
asteroid (Kaasalainen et al. 2001). However, in our relatively low S/N regime this method is not so useful. 
By binning our data over entire observational epochs, we average over the periodic flux variations due to 
the asphericity of the asteroid and solve for an ‘effective diameter’. Our interquartile mean provides us with 
statistically robust data at each viewing geometry. As a test, we did perform fits for 2 of our asteroids using 
all the observations separately, without binning, and the results were found to agree well with the results 
found using our binning process. For more on the assumption of sphericity in our model, see Section 

In our modeling we used Keplerian orbital parameters from the MFC, and found the absolute magni¬ 
tudes of the asteroids using the JPL Horizons web in terfac^ which were assigned an uncertainty of 0.3 
magnitudes, as was done in Mainzer et al. (2011a) and Mainzer et al. (2011cI. 


3. Methods 


Markov chain Monte Carlo (MCMC) methods sample probability distributions by constructing a Markov 
chain in state space which converges on the desired equilibrium distribution. A discussion of the mathematics 
behind the algorithm is beyond the scope of this paper (see jMackay 2003), but the method is often used in 
astronomy to sample posterior probability distributions of free parameters in a model given some data. It is 
useful to think of a Markov chain as a biased random walk, where the bias is such that the ‘walker’s’ steps 
converge to the desired probability distribution. Here, this is accomplished by defining a likelihood function 
using the familiar statistic which has as its equilibrium distribution the likelihood of a given parameter 
vector H being the ‘true’ parameter vector. We define: 




L[S] = Ke-sx [“1 

I fdata,i fmodel, 
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( 2 ) 

( 3 ) 


where i indexes the data points, each of which has a flux fdata,i, an uncertainty on that flux and a 
time of observation ti. k is a normalization constant which may be ignored for our purposes since MCMC 
methods evaluate only L[Hi]/L[H 2 ] to determine the acceptance or rejection of the next parameter vector. 
We have assumed a diagonal covariance matrix on the data in our equation for simplicity, which should be 


a good approximation. We used the emcee package for our MCMC analysis (Foreman-Mackey et al. 2013). 


emcee provides an ‘ensemble sampler’ which is afHne-invariant and utilizes a large number of ‘walkers’ to 
efficiently explore and sample parameter space, while employing parallelization to reduce the computational 
time needed for sampling. Affine-invariance ensures that the performance of our MCMC is not affected 


by correlations between our parameters causing anisotropic probability distributions (Goodman & Weare 


2010 ). 


Our thermophysical model has five free parameters: 


^Pseudocode: 
do M times 

take N samples from a unit normal distribution 

calculate the standard deviation of the combined second and third quartiles 
calculate the average standard deviation 

divide by the standard deviation of N/2 samples from a unit normal distribution of 1 / y/ N/2 
This calculation produces C ~ 0.77, indicating that the data in the first and fourth quadrants which are thrown out are 
nonetheless adding information to our statistic and so need to be accounted for. 

^http: / / ssd.jpl.nasa.gov/?horizons 
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• ip - The RA of the spin axis of the asteroid 

• 9 - The Dec of the spin axis of the asteroid 


01 - The dimensionless thermal parameter of Spencer et al. (1989) computed at a distance of 1 AU 


• e - The emissivity of the asteroid surface 

• D - The effective spherical diameter of the asteroid 

The thermal parameter is a measure of the importance of thermal inertia on the temperature. It is 
defined as 0 = , where k is the thermal conductivity of the regolith, p is the density, C is the specific 

heat per unit mass, is the rotational frequency of the asteroid, a is the Stephan-Boltzmann constant, 
and To is the equilibrium temperature of a fac et on the asteroid oriented toward the Sun. The combination 
r = y/npC is known as the thermal inertia (Winter & Krupp 1971). The equilibrium temperature is 

To = ) where A is the Bond albedo defined below, Lq is the solar luminosity, and r is the 

distance between the asteroid and the Sun. The parameter we fit is the thermal parameter when the asteroid 
is at a distance of 1 AU. For each data point we transform 0i as 0(r) = (r/lAU)^/^0i, using the value 
of r during that observation, in order to get the thermal parameter that should be used in calculating our 
model fluxes for that data point. We approximate the emissivity as independent of wavelength. Since we 
approximate the asteroid as spherical, the diameter we fit is an effective diameter. 

At each MCMC step, we use the diameter and the absolute magnitude of the asteroid to calculate a 
Bond albedo A for use in our calculations. The relation used is: 


A = q 


1329 km 10-^/^ 
D 


o r I {a) ■ 
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where q is known as the ‘phase integral’, and accounts for the fact that the brightness I of the asteroid 
may depend upon the phase a. q may be evaluated analytically using a full model of the reflectivity of an 
object’s surface, or an empirical number may be calculated. We use the q = 0.384 value used in |Mainzer et| 
al. (2011b), which is based on the lAU Commission 20 recommended value of the slope parameter for the 


magnitude-phase relationship of 0.15 in the H — G model of Bowell et al. (1989). Since the true value of q 


is not well known, this could be a possible cause behind cases where the model fails. An incorrect value of 
q or H could affect the probability distribution for the diameter in a complicated way, since D enters our 
model both by itself and through A. Where we have data at all four wavelengths, we fit for the albedo in 
the two shortest wavelengths, making the simplifying assumption that the albedo at 3.4/rm is the same as 


the albedo at 4.6^m (Mainzer et al. 2011bI. 


We used priors on the angular elements to enforce the modular definition of the spin axis, and used 


a uniform prior of 0.9 < e < 1.0 to force a physically sensible solution (Salisbury et al. 
emissivity is in general poorly constrained. 


1991), since the 


Wright ( 2007[ ) gives a de- 


The computationally-heavy step of our MCMC is computing fmodel, 
scription of the SCREAM algorithm. Here, we present details of the specific numerical methods used in 
our pilot study. As a general note, a large amount of time can be saved by taking full advantage of array 
operations and linear algebra in the following calculations and manipulations. 

Our primary concern is to find the temperatures of facets of representative craters at different latitudes 
on the surface of the asteroid over time by solving the power balance equation. For this we must calculate 


the sources and sinks of power for each facet of each crater. As Wright (2007) lays out, we account for the 


incident solar flux, the blackbody radiation of the facet, solar reflection and blackbody radiation from other 
visible facets, and vertical heat conduction into and out of the surface of the facet. 

First, from the Keplerian orbital elements of the asteroid and the time of the observation, the positions 
of the Sun and Earth are calculated in the asteroid-centric frame. We then go about creating the facets of 
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our representative crater at a particular latitude and local time. We utilize a cylindrical projection to go 
from facets in the 2d Cartesian square {x,y) G [—1,1] x [—1,1] to the 3d spherical cap geometry we require 
via the transformation: 


Or, 


arcsin(y) + — 


ip = TTX 


where the radial size is of no consequence since these are only representative craters, and Omax = 7r/4 was 
used in our modeling. This transformation maps lines parallel to the x axis to azimuthal circles around the 
spherical cap, and lines parallel to the y axis to great circles through the zero polar angle point. It also 
makes it simple to ensure our facets have equal surface area through the correct choice of bounding boxes 
before projecting, which we can calculate by changing variables under our map. 

Chains with Omax as an additional free parameter were run for several asteroids. Minimum went down 
only marginally in each case, and the posterior probability distributions were found to be wide with a peak 
near 7r/4. Thus Omax = 7i'/4 was used to model the generalized surface roughness which causes the beaming 
effect. Other works have taken similar strategies of removing Omax as a free parameter, with Hansen (19771 
assuming Omax = tt/S and Rozitis & Green (20111 performing fits using a few different angles, for example. 

We must choose how fine our discretization of the spherical cap will be. We decide for simplicity’s sake 
to use the same number n of facets in both the polar and azimuthal directions, for a total of facets. Fully 
optimized, our code is ~ O(n^), so we have a trade-off between computational expense and model accuracy. 

Slightly different temperature distributions are found for models with greater numbers of facets due to 
partially-shadowed facets on crater rims being broken into some shadowed facets and some lit facets. While 
the fraction of the facets which are lit is the same in each model, a model with more facets more accurately 
calculates the average surface normal vector of lit facets in the crater. Since this causes temperatures to 
change, it creates a linear dependence on frequency of the radiated flux change due to the exponential factor 
in Planck’s law. We note that the fluxes do approach some limit as the facet number increases, and 32^ 
facets gets close. If we go down to a model with 12^ facets, the cumulative loss in accuracy is usually ~ 6% 
at the high frequency end, which is well within the minimum 0.1 In 10/2.5% ~ 9% uncertainty we have 
provided for our data points. Since this is a proof-of-concept study and our computational resources were 
limited, we performed our calculations with 144 facets. 

It’s important to note that in general the error from using fewer facets depends upon the phase angle 
of the observation and the parameters used in the modeling. Moreover, in some cases, this error can be 
higher than the typical < 9% uncertainty. After we found best-fit parameters for each of our asteroids, we 
rechecked the loss in accuracy with these parameters for the specific viewing geometry of each data point. 
In each case, the flux errors were less than or comparable to the uncertainty we added to our data points, 
and corresponded to changes in predicted brightness temperature for the asteroid of less than IK for all 
frequencies. This level of accuracy was deemed sufficient for this proof-of-concept study. 

Once we have chosen our facet division, we set up the surface area integral and set it equal to one facet 
area: 


A = 


+ i 27r{l- COS Omax) 

/ r sm 0 dO dtp = — 
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144 


We now change variables under the aforementioned map taking r = 1 for ease of evaluation, and we 
require {xi,yj}j^ satisfy the following condition: 


(1 COS^^aa,) 

144 


= (a;i+i - Xi) 

(arcsin(j/j) + arcsin(j/j_|_i) -|- tt) 
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Solar Visibility of Crater Facets 



Figure 1: A visualization of the solar visibility of some representative craters over the face of a spherical 
asteroid. The lighter parts have a direct line of sight to the Sun, and the dark parts are in shadow. The 
sub-solar point is marked with a star. 


We take Xi = for azimuthal symmetry, and then can easily solve for the {yj} using computational 
methods. 

Defining the facets in this manner places a spherical cap at the origin opening downward, so we must 
maneuver it into the correct position. We first place it on the midnight longitude at the point on the surface 
of the asteroid directly opposite from the Sun, and then rotate it to the equator around the cross product 
of the spin axis and the location of the Sun. Then we can simply rotate it to the right latitude, and then 
around the spin axis to the correct local time. 

With the facets in place, we can compute the cosines of the angles between the surface normals and the 
Earth and Sun. We can now also compute the ‘solar visibility fraction’ of each facet - that is, the fraction 
of the facet which has a direct line of sight to the Sun, and thus incident solar flux. See Figure for an idea 
of the geometry involved. 

This is accomplished by considering a finer breakdown of each facet into 100 ‘pixels’, and calculating the 
point of intersection between the sphere of the crater and the line through each pixel and the Sun. For a 
sphere defined as |x — c|^ = 1^ with x the points on the sphere, c the center, and the radius taken to be 1; 
and a line defined as x = o -F with o the pixel location, i the unit vector between the pixel and the Sun, 
and s a free variable over the reals, we can calculate the distance s from the pixel at which the intersection 
lies as s = 2.£ • (c — o). If the point given by o -F s£ is farther from the Sun than the pixel in question (if 
s < 0) or the point is closer to the bottom of the crater than -y/l -F sin^ ^max, the maximum extension for 
our spherical cap, then the ray between the pixel and the Sun does intersect the crater, so the pixel does 
not see the Sun. If not, then there will be solar flux incident upon that pixel. We do this for all of the 100 
pixels in each facet, and then can calculate the total solar power incident upon the facet. 

Once we have calculated the solar visibility fractions of each facet of a crater at a given latitude at each 
time step as it rotates around the asteroid, we must determine the vertical heat conduction before we can 
solve the power balance equation. 

The vertical heat conduction coefficients G couple the temperature of a facet at one time to the tem¬ 
perature of the same facet at different times. We calculate them by computing the sum given in [Wright] 
(2007) Equation 14. Due to the dependence here, we can cut off the sum after 100 terms with little 
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loss (< 1.2 %), which speeds up this step a fair amount. Once this step is complete, we have precomputed 
the coupling coefficients between the temperatures of all the facets of a crater at all times. We then need to 
solve this system of equations simultaneously for the temperatures of the facets. 

To do so, we create a function that takes in 32 x 144 = 4608 temperatures of crater facets at a given 
latitude over time, and returns the residuals of the power balance equation calculated by subtracting the 
RHS of Equation 21 in Wright (20071 from the LHS. Even though the equations are coupled and non¬ 
linear, they can still be solved using iterative numerical methods. The best method found was a Newton- 
Krylov solver implemented in the SciPy Python packag^ which uses the inverse Jacobian to iteratively 
minimize the residuals. This solver generally returns sensible solutions, and converges relatively quickly to 
the max(residuals) < 10“^ level in ~ 100 evaluations. Examples of the temperature distribution returned 
by the solver can be found in Wright (2007). 

After running this for each latitude, we can use the temperatures of each facet of each crater over the 
asteroid to calculate the total blackbody radiation emitted by considering each as an independent blackbody. 
We add to this the reflected flux from the Sun by multiplying the solar flux by the appropriate geometric 
factors we found earlier. 

While we can evaluate this predicted flux at any frequency, rather than integrating over WISE’s passband 
spectral responses ourselves we use the linear quadrature formulae given in Wright (2013). We compute the 
predicted flux at 2 or 3 frequencies per passband, and by summing these with the appropriate weights we 
obtain analogues for WISE flux observations which are accurate to within half a percent. 

Now we can calculate the statistic by comparing these predictions to the data, and then feed that 
to the MCMC. Theoretically, we could start the Markov chain Monte Carlo at any point in parameter 
space, let it run, and - once we’ve waited long enough - return to hnd that the output of the chain has 
converged exactly to the posterior probability distributions we want. However, between the large number of 
computations needed to evaluate the model and the dimensionality of our parameter space, this approach 
is computationally unfeasible. We instead assist emcee in finding the global minimum, and then use it to 
map out that minimum and give us the posterior probability distribution. 

We begin by finding a reasonable range for the asteroid diameter by calculating the with various 
diameter guesses and by visual comparison with the overall magnitude of the flux. This process could 
easily be automated by testing some set values of the other four parameters with a set of diameter values, 
and is only done to save computational time. Then we coarsely sample parameter space with emcee by 
running a chain beginning at a widely-distributed ‘sample ball’, which starts walkers covering essentially the 
entire parameter space in (/?, 0,e and 0i, as well as the conceivable range in D. After we are satisfied that 
the parameter space has been well-investigated and decent wells have been found, we restart the chain 
centered on these different wells with smaller walker balls to find the minima of the wells. After investigating 
any promising wells, we are reasonably satished that the deepest of them is the global minimum. Eor our 
final sampling, we start two chains at this minimum of 50 steps with 50 walkers each, and then throw away 
the first 10 steps from each chain to account for ‘burn-in’. 


4. Results 

Characterizations of the posterior probability distributions (PPDs) of the parameters for each of the 
asteroids can be found in Table [2] Estimates of the thermal inertias of the asteroids can be found in Table 

Only fits where the best per degree of freedom was less than 1.5 are included in these tables. The 
number of degrees of freedom is the number of observational epochs multiplied by the number of operational 
passbands in each observational epoch, less the five fit parameters in our model. Ten out of the twenty 
asteroid fits which were attempted yielded such fits. 

For some examples of good fits, the best fit solutions for the asteroids 2014 HJ129 and 2102 Tantalus are 
plotted in Figures and respectively. The PPDs returned by the MCMC can be seen in Figures and 


®http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.newton_krylov.html 
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For an example of a failed fit, the best fit solution found for the asteroid 1627 Ivar is plotted in Figurej^ 
While it is apparent from the figure and the best per degree of freedom found of 90.76 that the SCREAM 
does not model this asteroid well, we are unable to say anything concrete about why the SCREAM failed 
in specific cases. All that is apparent is that no set of asteroid geometric and thermal properties could be 
found which produced the correct flux curves that 1627 Ivar was observed to have. However, this asteroid 


is one of those Hanus (20151 found a good fit for using their varied shape thermophysical model, indicating 


that it is likely that the assumption of sphericity is behind the failure of the SCREAM in this case. 


5. Discussion 


Since the model makes a number of assumptions, it is not surprising that good fits were not found for all 
asteroids. No beaming effect takes place with a sufficiently smooth asteroid, so a less-cratered model would 
be more appropriate in that case, and could perhaps be addressed by including a crater-covered fraction 
parameter in our model. For slowly-rotating asteroids, the data which were binned in each observational 
epoch may not have covered an entire period, or may not have been uniformly distributed over the asteroid’s 
rotation. This could cause our flux interquartile means to poorly represent the effective spherical diameter. 
Asteroids which are very aspherical could have a similar problem. We note that the SCREAM can be used to 


calculate fluxes for asteroids of arbitrary convex shapes (Wright 20071. While this generalization is useful 


convexity is not always a good assumption, especially in the case of bifurcated asteroids (near 10% of radar- 
imaged near Earth asteroids larger than ~ 200m are bifurcated candidate contact binaries (Benner et al. 


2006)) and exotic shapes like that of 25143 Itokawa (Demura et al. 2006). A new method would be needed 
to model more complicated shapes. Attempting fits of aspherical convex shapes to asteroids for which the 
spherical model failed should be done in follow-up work, as should an extension of this model to include 
a parameter like a crater-covered fraction. The use of the SCREAM in the varied shape thermophysical 


model method of Hanus (20151 when optical photometric data are available may also provide a fruitful way 


of incorporating shape effects. 

After obtaining our results, we examined the characteristics of our asteroids as well as our data for trends 
which might explain for which asteroids good fits could not be found. No patterns in orbital characteristics 
were found. No good fits could be found for any of the four asteroids with absolute magnitude H brighter 
than 15.0. Additionally, the SCREAM failed for four of the five asteroids with the largest normalized 
lightcurve amplitudes Overall the average normalized lightcurve amplitude of the asteroids for which 
good fits were not found was 3.2, while for the asteroids for which good fits were found the average was 2.8. 
Together, these suggest that shape effects are the largest problem for the SCREAM, which is unsurprising. 
Note that the lightcurve amplitude may be underestimated in the case of slowly-rotating asteroids, so it’s 
possible that some of the other failures are due to asteroids with significant shape effects which are rotating 
slowly enough for our method of lightcurve amplitude evaluation to be unable to discern. 

We also compared our modeled diameters and albedos to fits of the NEATM to our asteroids which have 


been previously published by the NEOWISE team in Mainzer et al. (2011c) and Mainzer et al. (2014). See 
Table 1^ for a comparison of the two. The two measures are generally in good agreement. 

Two of the asteroids we fit have classified taxonomic types: 2102 Tantalus is Q type and (89355) 2001 
VS78 is Sr type, placing both in the roughly stony category (Neese, 2010). The relatively high albedos found 


for both of these asteroids are consistent with this classification. 

Having a further empirical method for evaluating the fit parameters would allow us to check the external 
consistency of the SCREAM and our characterization methods. Unfortunately, radar data do not exist for 
any of the asteroids we characterizecQ A follow up comparison once such data are available would be prudent. 


®Normalized lightcurve amplitudes were calculated for each observing epoch and each band by sorting the data points, 
subtracting the 3N/4th value from the N/4th value - where N is the number of observations in that epoch - and dividing by the 
interquartile-averaged uncertainty of that epoch. We then averaged this measure over all of the observing epochs and bands 
for each asteroid. 

^http: / / echo.jpl.nasa.gov/asteroids / index.html 


























6. Conclusion 


The SCREAM can be used to derive asteroid effective spherical diameter, emissivity, albedo, thermal 
parameter, and spin axis when multi-epoch thermal IR data are available. Future follow-up work is needed 
to determine the accuracy of our fits. 

The information this gives us about the asteroids’ regoliths is limited because the thermal parameter 
depends on both the thermal inertia and the rotational period. This is reflected in the large uncertainties 
in TableEven with the uncertainties, however, some broad inferences can be made about which asteroids 
have surfaces which might be covered in rocky debris (thermal inertia ~ 50 J/m^/K/s^/^) or sand 400 
J/m^/K/s^/^) versus which are bare rock (> 2500 J/m^/K/s^^^). Additionally, the Yarkovsky effect depends 
solely on the thermal parameter, so this characterization may help in learning more about the typical 
magnitude of the Yarkovsky effect, and thus its importance ( |Delbo’ et al. 2007). 

This modeling may be most useful in the data it gives us on asteroid spin axes. Only a few hundred spin 
axes have been determined, as the commonly-used thermal models do not solve for them, and those that 


have been determined are biased toward larger asteroids (Kryszczyiiska et al. 2007). While estimates have 


recently been made for the longitudes of many spin axes using Lowell Observatory lightcurves, this method 
requires asteroids which exhibit relatively large absolute brightness variation, and does not provide the 
asteroid’s obliquity (Bowell et al. 2014). Characterization of asteroid obliquities and latitudes for smaller 
asteroids will help confirm spin axis distribution anisotropies and determine the relative importances of 
collisional, tidal, and thermal processes in altering spin directions, which is of some theoretical interest 
(Vokrouhlicky et al. 2003). Before using our results for these purposes it would be prudent to do follow¬ 


up testing of our results, including using more computational power and testing models which incorporate 
additional parameters such as non-spherical shapes. 

As NEOWISE continues its three-year mission, the number of asteroids observed at multiple epochs will 
increase dramatically, and MCMC modeling using the SCREAM provides a method for characterizing these 
asteroids. 
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A summary of the data used in our fits. The MJD given is the average modified Julian date of data points which were mid-averaged at that 
epoch. The RA and Dec are the RA and Dec of the asteroid as seen from Earth, r and A are the distances between the asteroid and the Sun 
and Earth, respectively. Since WISE orbits the Earth, to the relevant accuracy A is also the distance between WISE and the asteroid. Fluxes 
and uncertainties in WISE’s four bands are all in Janskys. A flux of 0 is a non-detection in that band. The average signal to noise ratio (SNR) 
is calculated by dividing the flux in each band by the uncertainty in that band, and averaging across that observation. 




Table 2: Parameter Posterior Probability Distribution Characterization 
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+la 

76° 

7° 

38.1° 

0.32 

0.963 

1663 

0.56 

2001 VS78 

Median 

51° 

-18° 

19.4° 

0.18 

0.930 

1608 

0.43 


— la 

33° 

-34° 

12.2° 

0.09 

0.907 

1579 

0.32 


+la 

20.° 

17° 

45.8° 

0.61 

0.980 

1794 

0.54 

1999 GJ4 

Median 

7° 

1° 

34.8° 

0.38 

0.946 

1732 

0.41 


— la 

354° 

-12° 

30.3° 

0.17 

0.916 

1673 

0.30 


Posterior probability distribution characterization of the free parameters in the SCREAM applied to the 
10 well-characterized asteroids (out of the 20 attempted characterizations) as well as the geometric visible 
albedo, a derived parameter. When deriving a posterior probability distribution for the albedo, we needed 
to account for the uncertainty in external values we use to calculate it along with the diameter. We therefore 
built a PPD using values of H drawn from a Gaussian of width 0.3 magnitudes centered at the value from 
JPL Horizons. The ±lcr values are found by sorting the PPD for a parameter and taking the 50 ± 34% 
values. 
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Table 3: Thermal Inertia Estimates 



Thermal Inertia 
(J/m^/K/s^/^) 

2014 HJ129 

13901^20 

2009 UX17 

24501^7° 

2102 Tantalus 

^70+230 

^'^-240 

2000 PJ5 

5404f9° 

1998 SB15 

lOOOl^^o 

2000 RJ34 


2010 NG3 

OOU 2io 

2003 LC5 

1570l«« 

2001 VS78 

10lt^« 

1999 GJ4 

9on+200 

ZZU_i3o 


Estimates of the thermal inertia of the fit asteroids. The thermal inertia depends on the thermal parameter, 
the albedo, the emissivity, and also the rotational frequency of the asteroid. We have PPDs of the first 
three of these from our fits, but not of the rotational frequency. We created PPDs for the thermal inertia 
by sampling the rotational frequency from a uniform distribution between (24 hours)and (2 hours)”^, 
from the findings of Pravec (20081. The thermal inertia values we give here are the median and ±lcr values 
of the derived PPDs. 


Table 4: Comparison of SCREAM and NEATM Models 



SCREAM 
Diameter (m) 

NEATM 
Diameter (m) 

SCREAM 

Albedo 

NEATM 

Albedo 

2014 HJ129 

5io;26 

519 ± 190 

u.uoy_Q Q]^Q 

0.031 ±0.027 

2009 UX17 

306;i? 

309 ± 10 

u.u^o o.oii 

0.042 ± 0.008 

2102 Tantalus 

1503^;^? 

1810 ±214 

o-342;o;“ 

0.214 ±0.084 

2000 PJ5 

757;2? 

923 ± 10 

0.308;°;“ 

0.227 ±0.033 

1998 SB15 

363ti 

337 ± 26 

f) nc7+0.022 

U.UO/o.oi5 

0.062 ±0.014 

2000 RJ34 

3893^^^ 

4330 ± 100 

U.UU/_Q 

0.067 ±0.011 

2010 NG3 

1149;^^ 

1520 ± 40 

n 1 70+0.063 
U.l ( O o.o44 

0.100 ±0.022 

2003 LC5 

1667;2? 

1678 ±10 

U.UOZ Q Q^3 

0.048 ± 0.048 

2001 VS78 

1608;^^ 

2000 ± 400 

U.+ZO 0.107 

0.240 ±0.131 

1999 GJ4 

1732;“ 

1641 ± 53 

U.^UU_o 102 

0.453 ±0.087 


A comparison between the modeled diameters and albedos from the SCREAM and those from the NEATM, 
which were found in Mainzer et al. (2011c I and Mainzer et al. (20141. The albedos are the geometric visible 
albedos. Note that best-fit absolute magnitude estimates from the JPL Horizons system may change over 
time as more data becomes available, so the values used in earlier works and in this paper may be different. 
This would in turn affect the relationship between the diameter and the albedo. 
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Figure 2: Best fit solution for the asteroid 2014 HJ129 with 3 DOF and = 0.60. Blackbody radiation 
is the dashed black line, reflected solar light is the dotted black line, and their sum is the solid grey line. 
Data points are plotted with error bars, which may be too small to be seen. The fluxes are not normalized 
for different heliocentric distances between epochs, but appear as observed. 
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Figure 3: Posterior Probability Distributions from a Markov chain Monte Carlo fit of the SCREAM 
applied to the asteroid 2014 HJ129. Visible on the plots are the 2.5% to the 97.5% bins, except for the 
emissivity in which the entire allowable range is shown. The solid vertical line denotes the median value, 
and the dashed vertical lines denote the Tier values, as reported in Table [2| 
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Figure 4: Best fit solution for the asteroid 2102 Tantalus with 7 DOF and = 2.84. Blackbody radiation 
is the dashed black line, reflected solar light is the dotted black line, and their sum is the solid grey line. 
Data points are plotted with error bars, which may be too small to be seen. The fluxes are not normalized 
for different heliocentric distances between epochs, but appear as observed. 
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Figure 5: Posterior Probability Distributions from a Markov chain Monte Carlo fit of the SCREAM 
applied to the asteroid 2102 Tantalus. Visible on the plots are the 2.5% to the 97.5% bins, except for the 
emissivity in which the entire allowable range is shown. The solid vertical line denotes the median value, 
and the dashed vertical lines denote the Tier values, as reported in Table [2| 
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Figure 6: Best fit solution of the failed fit for the asteroid 1627 Ivar with 1 DOF and = 90.76. 
Blackbody radiation is the dashed black line, reflected solar light is the dotted black line, and their sum 
is the solid grey line. Data points are plotted with error bars, which may be too small to be seen. The 
fluxes are not normalized for different heliocentric distances between epochs, but appear as observed. 
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